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PREFACE 


The  work  described  in  this  report  was  conducted  by  J.  E. 
Edinger  Associates  for  the  U.  S.  Army  Engineer  District,  Savannah 
(SAS) ,  in  partial  fulfillment  of  Contract  No.  DACW21-80-C-0042. 
The  purpose  of  the  study  was  to  investigate  the  feasibility  of 
applying  the  Laterally  Averaged  Reservoir  Model  (LARM)  with 
estuarine  boundary  conditions  to  narrow,  stratified  estuaries. 

The  feasibility  study  was  jointly  sponsored  by  the  U.  S. 

Army  Engineer  Waterways  Experiment  Station  (WES)  and  the  SAS. 

The  WES  participation  was  sponsored  by  the  Office,  Chief  of  Engi¬ 
neers,  U.  S.  Army,  under  the  Environmental  Impact  Research 
Program  (EIRP) . 

The  study  was  managed  by  Mr.  Robert  Gladden,  SAS.  Mr. 

Ross  Hall,  WES,  was  the  EIRP  Principal  Investigator.  Dr. 

Billy  H.  Johnson,  Hydraulics  Laboratory,  WES,  provided  the  tech¬ 
nical  review.  General  supervision  was  provided  by  Mr.  Donald  L. 
Robey,  Chief,  Ecosystem  Research  and  Simulation  Division,  WES, 
and  Dr.  John  Harrison,  Chief,  Environmental  Laboratory,  WES. 

Mr.  John  Bushman  was  Technical  Monitor  for  the  Office,  Chief  of 
Engineers. 

Commanders  and  Directors  of  WES  during  the  preparation  of 
this  report  were  COL  Nelson  P.  Conover,  CE,  and  COL  Tilford  C. 
Creel,  CE.  Technical  Director  was  Mr.  Fred  R.  Brown. 

This  report  should  be  cited  as  follows: 

Edinger,  J.  E.  and  Buchak,  E.  M.  1981.  "Estuarine 
Laterally  Averaged  Numerical  Dynamics;  The  Develop¬ 
ment  and  Testing  of  Estuarine  Boundary  Conditions 
in  the  LARM  Code,"  Miscellaneous  Paper  EL-81-9, 
prepared  by  J.  E.  Edinger  Associates,  Inc.,  for  the 
U.  S.  Army  Engineer  Waterways  Experiment  Station, 

CE,  Vicksburg,  Miss. 
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CONVERSION  FACTORS,  METRIC  (SI)  TO  U.  S.  CUSTOMARY 
UNITS  OF  MEASUREMENT* 


Metric  (SI)  units  of  measurement  used  in  this  report  can  be 
converted  to  U.  S.  customary  units  as  follows: 


Multiply 

centimetres 

kilometres 

metres 


_BX_ 


0.3937 

0.5396 

3.281 


To  Obtain 


inches 

miles  (U.  S.  nautical) 
feet 


*  U.  S.  nautical  miles  (nm)  are  used  in  this  report  to 
describe  the  study  area,  define  model  segment  boundaries, 
and  designate  the  locations  where  data  were  collected  by 
the  Chesapeake  Bay  Institute. 
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ESTUARINE  LATERALLY  AVERAGED  NUMERICAL  DYNAMICS  : 


THE  DEVELOPMENT  AND  TESTING  OF  ESTUARINE 


PART  I:  INTRODUCTION 


This  report  presents  the  results  of  a  feasibility  study  of  the 
application  of  the  laterally  averaged  reservoir  model  (LARM)  with 
estuarine  boundary  conditions  to  the  Potomac  River  estuary.  The 
feasibility  study  was  jointly  sponsored  by  the  U.  S.  Army 
Engineer  Waterways  Experiment  Station  (WES)  and  the  U.  S. 

Army  Engineer  District,  Savannah. 

The  LARM  code  was  selected  for  the  following  reasons: 

a.  It  has  been  applied  successfully  to  a  number  of  lab¬ 
oratory  and  reservoir  cases  over  the  past  three  years 
both  within  and  outside  the  Corps. 

b.  The  code  is  quite  general  for  the  hydrodynamic  equa¬ 
tions  and  can  handle  features  such  as  layer  additions 
and  subtractions  and  irregular  bottom  topography. 

c.  The  formulation  is  spatially  implicit  and  allows  for 
efficiently  long  computational  time  steps  in  deep 
waterbodies . 

d.  The  code  is  presently  being  extended  to  compute  the 
transport  of  many  interacting  water  quality  constit¬ 
uents  through  the  water  quality  transport  module  ( WQTM) . 


Although  the  LARM  code  applies  to  reservoirs,  the  original  basis 
of  the  formulations  was  numerical  estuarine  dynamics  (Edinger,  1974; 
Edinger  and  Buchak,  1975). 


The  Potomac  River  estuary  was  selected  for  the  feasibility 
study  because  of  the  availability  of  pertinent  data  and  exist¬ 
ence  of  comparative  numerical  model  studies.  The  Potomac  is  a 
classical  coastal  plain  estuary  exhibiting  salinity  stratifica¬ 
tion  and  baroclinic  circulation  in  the  lower  reaches.  It  has 
detailed  geometry  data  readily  available  of  the  type  required 
by  LARM  (Cronin  and  Pritchard,  1975).  It  has  been  studied  in 
conjunction  with  various  tidal  and  hydrodynamic  models  (Rives, 

1973;  Blumberg,  1975;  Elliott,  1976;  Wang  and  Kravitz,  1980). 

There  is  a  short-term  extensive  data  base  of  horizontal  tidal  veloc 
ities  and  salinity  collected  especially  for  comparison  to  models. 
The  Potomac  is  similar  to  the  James  River  estuary  where  Pritchard 
(1956,  1967)  evaluated  the  distribution  of  vertical  velocity 
and  vertical  mixing  that  can  also  be  compared  to  model  results. 

The  Potomac  presents  an  ideal  situation  for  estuarine  model 
testing  and  demonstration. 

Development  of  LARM  for  estuaries  began  with  a  re-examination 
of  the  implicit  numerical  formulation  of  laterally  averaged  water- 
body  dynamics  and  the  characteristics  of  the  LARM  code.  An 
investigation  was  made  of  the  dynamics  at  the  downestuary  boundary 
where  the  velocity  profile  was  computed  using  the  specified  tide 
height  and  salinity  distribution.  The  only  terms  relying  on 
auxiliary  relationships  were  the  vertical  and  horizontal  dispersion 
processes.  The  vertical  mixing  processes,  in  particular,  are  quite 
important  in  narrow,  deep,  and  stratified  estuaries;  hence  an  exten 
sive  investigation  was  made  of  their  properties  for  use  with  LARM. 
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The  model  set  up  and  testing  included  the  steps  of  estab¬ 
lishing  the  model  geometry,  choosing  the  numerical  boundary  and 
initial  conditions,  validation,  and  initial  simulations  and 
results.  The  validation  step  tested  the  code  with  an  estuarine 
boundary  for  simple  cases  to  assure  that  there  are  no  computa¬ 
tional  anomalies.  The  results  of  the  initial  model  simulations 
were  summarized  and  examined  using  vector  plots  of  the  circulation 
in  the  estuary  throughout  a  tidal  cycle. 

The  model  results  were  compared  to  field  data  for  tide  height 
range  and  phase  along  the  estuary  using  the  data  in  Rives  (1973). 

The  velocity  over  a  number  of  tidal  cycles  at  different  depths 
and  velocity  profiles  at  different  times  in  the  tidal  cycle  were 
compared  to  the  observations  of  Elliott  and  Hendrix  (1976).  The 
vertical  velocity  and  dispersion  coefficient  variations  found 
with  the  model  were  compared  to  the  results  of  Pritchard  (1967) . 
These  comparisons  of  model  results  and  field  data  constitute  the 
model  verification.  The  previous  verifications  of  LARM  in  res¬ 
ervoirs  used  only  constituent  (temperature)  distributions  because 
few  field  measurements  existed  of  the  very  low  velocities  found 
in  these  waterbodies.  The  Potomac  estuary  comparisons  are  the 
first  velocity  verifications  of  the  LARM  code. 

Numerous  features  of  the  model  that  are  not  illustrated  by 
the  real-time  simulations  were  examined  using  sensitivity  analyses. 
These  included  the  rate  at  which  salinity  becomes  distributed 
throughout  the  estuary  when  initialized  from  zero  concentration, 
the  sensitivity  of  the  velocity  and  salinity  profiles  to  the 


vertical  dispersion  coefficient,  and  the  effects  of  wind.  The 
latter  was  examined  in  detail  because  wind  can  be  an  important 
force  driving  estuarine  circulation  in  real  situations. 

Along  with  comparing  model  results  and  field  data  to  demon¬ 
strate  model  verification,  the  report  elaborates  on  the  important 
features  of  estuaries  the  model  describes.  The  modeling  task 
should  eventually  be  simplified  so’  that  the  major  portion  of 
the  effort  can  be  dedicated  to  assembling  and  examining  the  time- 
varying  boundary  data  and  other  input  data  including  tide  record 
analysis;  boundary  salinities;  and  wind  speed,  direction,  and 
wave  height  estimates  for  dispersion. 

The  code  is  recommended  for  those  cases  where  longitudinal 
and  vertical  gradients  occur  and  where  lateral  homogeneity  can 
be  assumed. 
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PART  II:  ESTUARINE  DYNAMICS  IN  LARM 


r 

! 


The  numerical  description  of  estuarine  dynamics  is  formu¬ 
lated  from  the  basic  equations  of  fluid  motion .  A  discussion  of 
the  formulation  for  the  laterally  averaged,  vertically  mixed, 
and  sectionally  homogenous  estuary  cases  is  given  in  Edinger 
and  Buchak  (1980).  Details  of  the  laterally  averaged  relation- 

* 

ships  are  presented  here  to  show  the  computational  basis  of  the 
LARM  code  and  its  application  to  estuaries. 

The  LARM  reservoir  applications  differ  from  estuary  appli¬ 
cations  only  in  the  specification  of  boundary  conditions.  The 
LARM  code  also  incorporates  many  computational  algorithms  such 
as  generalized  bottom  coordinates,  addition  and  subtraction  of 
top  layers,  upstream  cell  addition,  and  input/output  routines  that 
broaden  its  use  in  application  to  estuaries. 

The  estuarine  boundary  condition  of  the  varying  tide  and 
salinity  at  the  estuary  mouth  requires  that  the  model  compute 
the  fluxes  into  and  out  of  the  estuary  from  numerical  dynamics. 
Certain  approximations  are  implied  in  the  numerical  boundary 
conditions  which  may  affect  results. 

Once  the  hydrodynamics  are  formulated,  there  are  only  two 
unknown  internal  parameters.  These  are  the  longitudinal  and 
vertical  turbulent  dispersion  parameters.  In  some  estuaries,  the 
latter  is  more  important  than  the  former  and  more  is  known 
about  its  formulation.  These  parameters  are  examined  to  de¬ 
termine  the  formulations  for  use  in  the  estuary  version  of  LARM. 
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Laterally  Averaged  Waterbody  Dynamics 


The  laterally  averaged  equations  of  fluid  motion  as  they 
apply  to  estuaries  can  be  derived  from  the  three-dimensional 
equations  of  fluid  motion  as  illustrated  in  Edinger  and 
Buchak  (1980).  There  are  six  unknowns  and  six  equations  in¬ 
cluding:  (1)  the  free  water  surface  elevation,  n;  (2)  the  pres¬ 

sure,  P;  (3)  the  horizontal  velocity,  U;  (4)  the  vertical  ve¬ 
locity,  W;  (5)  the  salinity,  S;  and  (6)  the  density,  p.  The 
six  equations  are:  (1)  the  free  surface  wave  equation;  (2)  the 
hydrostatic  pressure;  (3)  horizontal  momentum;  (4)  continuity; 
(5)  constituent  transport;  and  (6)  an  equation  of  state  relating 
density  and  salinity. 

The  laterally  averaged  equations  of  fluid  motion  and  trans¬ 
port  are  the  horizontal  momentum  balance: 


9  UB  ,  9  UUB  .  9WUB  1  9  BP  ,  3  (BA  3U/3x)  .  B3tZ 

~Tt  +  T3T  +  =  ”  p  S3T  +  Tk  x  ~Jz~ 


(1) 


where  B  is  the  estuary  width  as  a  function  of  x  and  z;  and  U  and 
W  are  the  laterally  averaged  horizontal  and  vertical  velocity 
components.  The  vertical  equation  of  motion  reduces  the  hydro¬ 
static  approximation: 


The  equation  of  continuity  becomes: 


(2) 


3  UB  +  |_WB  _  qB 


3  x  3  z 

where  q  is  the  side  or  tributary  inflow  per  AxAzB  volume.  Ver¬ 
tically  integrated  continuity  gives  the  free  water  surface 
relationship  of 


(3) 
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(•h 

3B  n  _  _3_  UBdz  - 

3t  3x 

n 

where  B  is  the  time  and  spatially  varing  surface  width  and  n  is 
the  free  water  surface  elevation.  The  constituent  transport  for 
salt  becomes: 

3BS  +  +  3SS  .  |_(BDxas/3x)  .  i .  sq  B  CS) 

where  Sq  is  the  tributary  source  flux  per  AxAzB  volume,  and 

P  =  R(S)  (6) 

is  the  equation  of  state.  Each  additional  constituent  such  as 
suspended  sediment,  temperature,  or  dissolved  oxygen  has  a 
balance  similar  to  Equation  5  with  specific  source  and  sink 
terms.  The  equation  of  state  is  similarly  modified  for  constit¬ 
uents  such  as  temperature  and  salinity  that  have  a  significant 
effect  on  density. 

Equations  1  to  6  constitute  six  equations  to  be  solved 
for  the  six  unknowns  of  U,  P,  W,  n,  S  and  p.  Lateral  averaging 
eliminates  the  lateral  momentum  balance,  the  lateral  velocity 
component,  and  the  Coriolis  acceleration.  The  computational 
problem  is  reduced  to  six  equations  in  six  unknowns  and,  most 
important,  to  two  coordinate  directions.  The  reduction  to  two 
coordinate  directions  is  the  main  feature  that  reduces  computa¬ 
tional  time  and  storage  over  the  three-dimensional  case. 

The  laterally  averaged  horizontal  pressure  gradient  in  the 
horizontal  momentum  balance  is  the  density  driving  force.  It 
can  be  expanded  to 


h 

qBdz  (4) 

n 
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(7) 


3 BP  =  B3P  +  P3B 
3x  3x  3x 

The  second  term,  P3B/3x,  represents  the  static  force  of  the 
fluid  on  the  x  projection  of  the  lateral  boundary  which  in  turn 
is  cancelled  by  the  force  of  the  boundary  on  the  fluid.  Thus, 
B3P/3x  represents  the  internal  fluid  horizontal  pressure  gradient. 
The  horizontal  pressure  gradient  is  evaluated  from  Equation  2 
to  give 

f  z 


gB  +  gB 


( 3p/3x)dz 


(8) 


n 


at  any  depth  z.  The  horizontal  pressure  gradient  is  divided 
into  the  two  components  of  the  surface  slope  and  the  vertical 
integral  of  the  horizontal  density  gradient.  The  first  term 
is  known  as  the  barotropic  gradient  and  the  second  as  the 
baroclinic  gradient.  The  horizontal  density  gradient  is  the 
major  driving  force  for  the  density  circulation  exhibited  in 
coastal  plain  estuaries. 

The  basic  characteristics  of  the  longitudinal  and  vertical 
free  water  surface  hydrodynamics  can  be  examined  through  evalu¬ 
ation  of  the  water  surface  relationship,  Equation  4.  The  ver¬ 
tical  integral  of  the  horizontal  flow  required  in  Equation  4 
can  be  determined  from  the  algebraic  forward  time  difference  of 
the  local  acceleration  of  horizontal  momentum  in  Equation  1. 
Formulation  of  the  forward  time  difference  of  UB  is  the  first 

step  in  evaluating  the  numerical  equations.  It  gives 

z 

U ' B '  =  UB  +  gBAt  3n / 3x  -  gBAt  (3p/3x)dz  +  F  At  (9) 

P  * 

n 
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where  Equation  8  has  been  substituted  for  the  horizontal 


pressure  gradient  and  Fx  is 


F  -  i- (BA  9U/8X)  -  K-DUB),  _  |TOE  Ifl 

x  9x  x  °x  3z  3x 


(10) 


The  vertical  integrals  of  the  various  terms  in  Equation  9  can 
be  further  evaluated  for  insertion  into  the  vertical  integral 
of  the  flow  required  in  the  free  water  surface  balance, 
Equation  4. 

The  vertical  integral  of  the  horizontal  pressure  gradient 
can  be  evaluated  from  Equation  8  to  give 
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The  first  term  on  the  right-hand  side  results  from  the  fact 
that  9n  /  3  x  is  a  function  only  of  x  and  is  constant  over  z. 

The  integral  of  width,  B,  over  depth  is  the  total  cross-sectional 
area  across  which  the  surface  slope  contribution  to  the  hori¬ 
zontal  pressure  gradients  acts.  The  second  term  is  the  force 
due  to  the  horizontal  density  gradient. 

The  vertical  integral  of  the  horizontal  shear  stress  can  be 
expanded  from  the  derivatives  of  3Bxx/3z  to  give 


B3- 


Tir  dz  =  Bh  Th 


B  T  - 
n  n 


3B  , 

Tx3i  dz 


(12) 


The  first  term  is  the  bottom  shear  evaluated  at  z=h  and  can  be 
evaluated  from  bottom  velocity  friction  relationships.  The  sur¬ 
face  shear,  Bn  xn  is  the  surface  wind  shear  component  parallel 
to  the  x  axis.  The  third  term  is  the  wall  or  bottom  shear  due 


12 
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to  the  horizontal  projection  of  the  sloping  sides  of  the  water- 
body  (3B/3z).  It  can  be  evaluated  as  bottom  shear  over  the 
projected  width  3B  at  each  elevation.  The  internal  velocity 
shear  cancels  out  of  the  vertical  integration. 

Collecting  the  various  terms  of  Equation  9  into 
Equation  4  gives  the  surface  elevation  equation  of 
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(13) 


With  the  n  or  surface  coordinate  terms  collected  on  the  left 
hand  side.  Equation  13  is  the  water  surface  equation  of  the 
integrated  waterbody.  Equation  13  is,  therefore,  a  numerical 
form  of  the  frictionally  dampened  long  wave  equation  for  an 
irregular  geometry,  stratified  waterbody. 

For  the  laterally  homogenous  estuary,  Equation  13  can  be 
simply  evaluated  implicitly  from 
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Where  n ' ^  is  the  new  time  level  value  of  the  surface  elevation 
at  a  finite  x  location  i;  is  the  total  cross-sectional  area 
between  i  and  i+1  locations  of  n '  .  and  n  ’  •  , ;  n  •  is  the  water 
surface  elevation  at  the  previous  time  step;  and  is  the  sum 
of  terms  on  the  right-hand  side  of  Equation  13  except  for 
lateral  inflows  and  outflows.  Equation  14  is  a  spatially 
implicit  surface  relationship  that  eliminates  the  gravity  wave 
stability  criterion.  It  can  be  evaluated  on  each  time  step  using 
the  efficient  Thomas  algorithm  for  a  tridiagonal  martrix.  Devel¬ 
opment  and  use  of  Equation  14  in  laterally  averaged  spatially 
implicit  hydrodynamics  have  been  presented  in  Hamilton  (1975) 
and  in  Edinger  and  Buchak  (1975,  1980). 

The  numerical  procedure  for  solving  for  the  six  unknowns 
on  each  time  step  is  to  compute  first  the  water  surface  elevations 
from  Equation  14  and  to  obtain  the  horizontal  velocity  compo¬ 
nents  from  Equation  9.  The  vertical  velocity  component  is 
found  from  continuity,  Equation  3,  and  the  salinity  distribution 
from  the  constituent  balance,  Equation  5.  The  water  surface 
elevation  equation  essentially  results  from  the  simultaneous 
algebraic  substitution  and  solution  of  horizontal  momentum, 
Equation  1,  and  vertically  integrated  continuity,  Equation  4, 
giving  U  and  n  simultaneously.  This  substitution  makes  the  solu¬ 
tion  spatially  implicit  in  q  and  U  at  the  same  time  level, 
through  Equation  14,  and  eliminates  the  Courant  gravity  wave 

speed  criterion  that  Ax/At  >  /gH  which  leads  to  short  corn- 

max 

putational  At  in  deep  waterbodies. 
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With  the  laterally  averaged  equations  of  motion  expressed 
in  an  algebraic  form,  it  is  necessary  to  devise  a  finite  dif¬ 
ference  coding  for  numerical  computations.  The  coding  not  only 
includes  the  finite  difference  form  of  the  equations  but  also 
the  logic  and  algorithms  needed  to  carry  out  the  computations. 
These  procedures  have  been  developed  in  the  LARM  code. 

The  Numerical  LARM  Code 

The  variables  are  located  on  a  physical  space  and  computa¬ 
tional  grid  as  shown  in  Figure  1.  It  is  called  a  space 
staggered  grid  since  certain  of  the  variables  are  at  one  location 
and  certain  at  another.  There  is  a  rational  basis  for  choosing 
the  grid  locations  which  can  be  seen  by  using  imaginary  control 
volumes  about  a  point. 

The  constituent  concentration  S  is  surrounded  by  a  cell  that 
has  the  U  and  W  at  the  boundaries.  Therefore,  the  U  and  W  can 
transport  S  into  and  out  of  this  cell  with  no  spatial  averaging 
to  determine  the  change  in  S  over  time.  Similarly,  the  W  is 
computed  for  the  same  volume  using  the  U's. 

The  velocity,  U,  is  surrounded  by  a  cell  with  the  water  sur¬ 
face  elevations,  n ,  and  densities  known  at  either  end.  U  is 
computed  from  horizontal  gradients  of  the  surface  slope  and  den¬ 
sity  with  no  spatial  averaging  of  the  primary  variables.  Ad- 
vection  of  momentum  into  and  out  of  the  cell  does  require  spatial 
averaging  to  determine  the  fluxes  at  the  ends  of  the  cell,  but 
the  variable  U  being  computed  remains  centered. 

The  spatial  indices  of  Figure  1  define  how  the  finite 
differences  and  integrals  in  Equations  14,  9,  3,  and  5 
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are  formed.  There  is  extensive  spatial  averaging  connected  with 
each  of  the  terms  in  these  relationships  to  account  for  the 
spatially  varying  geometry.  This  averaging  is  detailed  in 
Edinger  and  Buchak  (1975,  1980)  and  need  not  be  expanded  on  here. 

The  free  water  surface  elevation,  n,  varies  within  the 
fixed  grid  in  Figure  1.  The  top  layer  thickness  varies  with 
n  and  the  water  surface  cannot  fall  below  the  computational 
grid  point.  When  +  *e  water  surface  is  rising,  it  can  in  prin¬ 
ciple  move  up  to  a r  elevation  and  the  top  layer  gets  deeper; 
the  computed  variables  would  simply  represent  lateral  averages 
over  the  greater  ^septh.  Letting  the  top  depth  increase  without 
refinement  leads  to  volume  errors  since  width  would  not  change, 
errors  on  surface  exchange  processes,  and  eventually  limited 
hydrodynamic  computations.  The  problem  is  severe  in  reservoirs 
where  water  surface  elevations  can  change  by  5  to  10  metres,*  yet 
vertical  resolution  is  required  to  1  metre.  It  is  particularly 
important  on  the  Savannah  River  which  has  a  mean  low  water  depth 
of  about  4  metres  and  tide  of  1  to  2  metres  with  a  spatial  reso¬ 
lution  required  of  0.5  metre.  This  is  the  top  box  problem. 

The  top  box  problem  has  been  solved  for  LARM  using  a  com¬ 
putational  algorithm.  The  procedure  is  that  on  a  rising  water 
surface,  where  the  elevation  in  all  cells  passes  a  specified 
level,  a  new  layer  is  added.  Since  the  variables  were  computed 
as  laterally  homogenous  before  the  new  layer,  they  will  then  be 

the  value  in  that  layer.  On  falling  water  surface  where  the 

*  A  table  of  factors  for  converting  metric  (St)  units  of 
measurement  to  U.  S.  customary  units  is  presented  on  page  3. 


elevation  in  all  cells  is  below  a  specified  level,  then  a  layer 
is  removed.  The  logic  to  carry  out  this  algorithm  is  quite  com¬ 
plex  and  is  incorporated  in  all  steps  of  the  computations. 

An  important  part  of  the  LARM  code  is  the  output  handling 
procedures.  The  basic  code  outputs  grids  of  the  main  variables, 
U,  W,  S,  n  at  specified  times.  For  reservoir  problems,  output 
of  results  once  per  day  was  often  satisfactory.  Estuaries  re¬ 
quire  at  least  hourly  output  and  more  frequently  if  tidal  average 
statistics  were  required.  One  part  of  the  study  has  been  to 
develop  procedures  for  retaining  results  and  the  routines  to 
perform  tidal  cycle  summaries.  It  can  become  very  costly  to 
store  all  the  output  results  for  any  period  of  time  without  the 
summary  routines. 

Estuarine  Boundary  Condition 

An  estuarine  boundary  condition  specifies  tide  height  (Z^) 
and  constituent  concentrations  (S.  .  ,  T.  )  at  the  ocean  bound- 

1  ,  K  X  ,  K 

ary,  and  computes  the  boundary  fluxes  (Ui_^  k).  The  basic  steps 
are  the  computation  of  the  upestuary  surface  elevations  using 
the  specified  tide  height  in  the  right-hand  boundary  side 
of  the  surface  wave  relation.  Equation  13,  and  then  computing 
the  boundary  U.  ,  ,  from  the  momentum  balance,  Equation  9, 

1“1  j  K 

across  the  i-1  to  i  interface. 

There  are  some  limitations  to  the  computation  of  the 
boundary  ^  k  from  the  space  staggered  computational  scheme 
across  the  boundary.  The  horizontal  momentum  computation  of 
Ui  k  in  general  depends  on  (1)  the  horizontal  pressure  gradient; 


(2)  horizontal  advection  of  momentum;  (3)  vertical  advection 
of  momentum;  (4)  horizontal  dispersion  of  momentum;  and  (5) 
vertical  dispersion  of  momentum  or  interfacial  shear  as  shown 
in  Equation  1.  Each  term  can  be  examined  separately  for  its 
evaluation  internally  and  at  the  boundary  to  show  the  limitations 
of  the  boundary  equations. 

The  horizontal  pressure  gradient  of  the  momentum  is  computed 
at  the  boundary  similar  to  the  interior  points.  For  the  space 
staggered  grid,  the  horizontal  pressure  gradient  internally 
falls  naturally  across  the  E\  ^  k  velocity  since  Z._^  and 
as  well  as  p.  ^  k  and  2  k  are  known  as  are  the  components 
of  the  pressure  computation.  The  same  terms  are  available  at 
the  estuary  boundary  since  Z.  and  p.  ,  are  specified,  and  Z  , 

1  1  ,  K  1  —  X 

and  p.  ,  ,  are  computed  internally. 

The  horizontal  and  vertical  advection  of  momentum  in  and 
out  of  the  box  require  having  the  advective  fluxes  with  which 
the  centered  U.  „  ,  is  transported,  and  also  having  the  concen- 

1“4j  ,  K 

tration  of  momentum  that  is  transported.  For  an  interior  lL_g  k> 
the  horizontal  transports  are  found  by  averaging.  Neglecting 
geometry  changes,  the  right-hand  transport  is  (IL  g  k  +  k)/2. 

This  transport  can  be  either  out  of  (+)  or  into  (-)  the  computa¬ 
tional  box.  When  out,  it  transports  U.  g  k  and  when  in,  it 
transports  U.  .  ,.  The  choice  of  the  quantity  being  transported 

1“1  (  K 

from  the  sign  or  direction  of  the  advective  fluxes  is  known  as 
upwind  differencing.  The  vertical  advection  of  momentum  is 
similarly  treated  with  the  vertical  fluxes  being  averaged  as 
(W.  „  .  +  W.  ,  )/2  to  be  centered  at  U.  0  .  and  the  momentum 
transported  chosen  by  the  sign  of  the  flux. 
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Advective  fluxes  of  momentum  are  not  readily  evaluated  at 


the  estuary  boundary.  There  are  no  U.  .to  the  right  of  Z. 

1  t  K  1 

unless  they  are  specified  from  data.  Similarly,  there  is  no 
W.  ,  since  there  are  no  U .  ,  for  its  determination  from 

1  ,  K  1  ,  K 

continuity . 

Horizontal  dispersion  of  momentum  is  computed  for  U.  0 

1“Z  ,  K 

from  the  horizontal  velocity  gradient.  For  an  interior  point, 
at  UA_2  k>  the  velocities  are  known  for  computation  of  the 
gradient.  For  ^  at  the  estuary  boundary  there  is  no 
known  U.  ,  for  the  gradient  computation.  Horizontal  dispersion 

1  ,  K 

of  momentum  is  not  readily  included  in  the  boundary  computation 

of  U.  . 

l-l ,  k 


Vertical  dispersion  of  momentum,  or  vertical  shear,  is  com¬ 
puted  from  the  vertical  gradient  of  U.  ,  using  the  vertical 

1  ,  K 

profiles  of  U.  .  centered  at  the  bottom  of  the  imaginary  com- 

1  ,  K 

putational  box.  It  can  be  computed  at  each  vertical  line  of  the 
velocity  components  including  those  of  the  boundary.  Bottom 
shear  due  to  the  changing  widths  in  each  layer  and  at  the  bottom 
are  also  readily  computed  at  the  boundary. 

At  the  estuarine  boundary,  the  horizontal  momentum  fluxes 
are  computed  from  the  horizontal  pressure  gradient,  the  internal 
vertical  shears,  and  boundary  friction.  Neglected  by  computation 
with  a  space  staggered  grid  are  the  horizontal  and  vertical  ad- 
vection  of  momentum  and  the  horizontal  dispersion  of  momentum. 

There  are  two  possible  effects  for  these  computational 
simplifications.  First,  the  momentum  terms  neglected  in  computing 
the  boundary  fluxes  could  become  balanced  out  in  the  next  upestuary 
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cell  which  then  becomes  the  real  boundary  cell.  Second,  since 
flux  at  the  boundary  is  balanced  almost  solely  by  internal  and 
bottom  friction,  the  solution  near  the  estuarine  boundary  could 
be  highly  dependent  on  the  vertical  momentum  dispersion  coeffi¬ 
cient  and  the  bottom  friction  factor.  These  coefficients  can 
be  treated  locally  for  the  boundary  cells  as  well  as  upestuary 
if  necessary. 

It  is  expected  that  there  will  be  some  artifical  circulation 
in  the  first  cell  upestuary  of  the  boundary  due  to  the  lack  of 
momentum  transport  by  advection  and  dispersion  across  it.  Com¬ 
parisons  of  observed  and  computed  velocities  upestuary  from  the 
boundary  will  indicate  if  this  circulation  has  any  major  effect 
on  the  flow  field. 

Vertical  Dispersion  Formulations 

Vertical  dispersion  of  momentum  enters  the  horizontal  momen¬ 
tum  balance,  Equation  1,  through  evaluation  of  the  horizontal 
shear,  t  .  The  shear  is  related  to  the  dispersion  process  as 

x  =  A  3U/3z.  Vertical  turbulent  transport  of  constituent 
z  z 

enters  the  constituent  balance,  Equation  5,  as  Dz3S/3z.  The 

terms  A  and  D  are  vertical  turbulent  dispersion  coefficients 
z  z 

for  momentum  and  constituent,  respectively,  and  are  the  only 
terms  in  laterally  averaged  numerical  relationships  that  require 
evaluat ion . 

The  source  of  the  turbulent  diffusion  terms  is  the  temporal 
averaging  that  is  assumed  to  apply  to  the  hydrodynamic  relationships 
to  give  time- smoothed  results  from  fluctuating  turbulent  flows. 
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The  arithmetic  of  this  step  is  to  represent  an  instantaneous 
value  of  velocity,  u,  by  a  short  time-averaged  mean  value  U  and 


the 

fluctuation  about 

the  mean. 

,  u ' ,  as 

u  = 

U  +  u' 

(15) 

and 

similarly  for  the 

vertical 

velocity  and  salinity 

w  = 

W  +  w' 

(16) 

s  = 

S  +  s' 

(17) 

Averaging  the  vertical  advection  of  momentum  term,  3UW/9z, 
and  the  vertical  advection  of  salinity,  9WS/ 9z  gives 

UW  =  U  W  +  <u'  w’>  (18) 


and 


WS  =  W  S  +  <w'  s'  >  (19) 

where  the  overbar  ( — )  signifies  a  time  average  of  the  product 
terms.  The  momentum  and  constituent  transport  balances  are  for 
the  mean  values,  U,  W,  S,  and  their  averaging  result  in  the 
average  cross  product  fluctuation  terms  <u'  w’>  and  <w'  s’> 
which  in  turn  represent  the  vertical  turbulent  fluxes  of  hori- 
zc ntal  momentum  and  salt. 

A  description  of  the  turbulence  processes  in  a  numerical 
model  is  complicated  by  the  relationship  between  the  model  integra¬ 
tion  time  step,  At,  and  the  time  period  of  averaging  implied  in 
Equation  18  and  19.  Presumably,  the  mean  values,  U,  W,  S, 
are  an  average  over  the  integration  time  step  and  the  fluctu¬ 
ations  about  the  mean,  u',  w' ,  and  s’,  are  measured  many  times 
within  this  time  step.  A  short  integration  time  step  of  a  few 
seconds  to  one  or  two  minutes  approaches  the  frequency  of  the 
fluctuations,  yet  the  model  is  not  actually  computing  turbulent 
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fluctuating  velocities.  Rather,  the  model  is  iterating  between 
the  time  limits  of  the  boundary  data  to  the  next  solution  point. 
One  advantage  of  an  implicit  solution  with  long  time  steps  is 
that  the  computed  U,  W,  and  S  at  the  end  of  each  time  step  rep¬ 
resents  a  mean  over  a  period  for  which  the  average  of  the 
fluctuations  apply. 

A  second  difficulty  of  describing  the  turbulent  transport 
processes  for  use  in  a  numerical  model  is  relating  them  to  the 
scale  of  the  grid.  Although  the  turbulent  transport  terms  are 
properly  multiplied  by  the  width  of  lateral  averaging,  the  mean 
velocity  and  salinity  in  the  gradient  terms,  3U/3z  and  3S/3z, 
apply  across  a  segment  area  BAx  where  B  is  of  the  order  of 
hundreds  of  metres  and  Ax  is  of  kilometres.  The  width  average 
is  implied  when  computing  U,  W,  and  S  but  not  the  length  average. 
Thus,  there  is  a  component  of  spatial  variation  contained  in 
<u'  w’>  and  <w'  s’>,  that  varies  over  the  computa*.  *e»ial  griu 
length,  Ax. 

The  average  cross  product  fluctuation  terms  which  make  up 
the  vertical  turbulent  transport  process  are  related  to  the  mean 
flow  and  concentration  fields  predicted  by  a  model  as 


<u'  w'>  =  -  Az  3U/3z 

(20) 

<w'  s’>  =  -  Dz  3S/3z 

(21) 

The  problem,  therefore,  reduces  to  evaluating  Az  and  Dz  in 
terms  of  the  computed  flow  and  density  fields. 

In  an  estuary  that  can  stratify,  the  vertical  dispersion 
coefficients  are  affected  by  the  level  of  turbulence  and  the 
degree  of  stratification.  They  can  be  represented  as 
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(22) 


A  =  A  F(Ri) 

Z  Zo 

where  A  is  a  function  describing  the  level  of  turbulence 
Zo 

related  to  velocity,  depth,  etc.,  and  Ri  is  the  Richardson  number 
related  to  stratification.  Almost  every  formulation  of  Az  that 
is  available  can  be  examined  in  the  form  of  Equation  22. 

The  Richardson  number  is  an  important  concept  in  describing 
vertical  mixing  in  a  stratified  waterbody.  It  is  defined  as 


Ri 


I  l£ 

p  3z 


2 


(23) 


which  is  the  ratio  of  the  buoyant  or  potential  energy  to  the 

kinetic  energy  being  dissipated  at  a  point  in  the  water  column. 

The  greater  the  stratification,  as  indicated  by  a  large  density 

gradient  3p/3z,  the  more  a  given  amount  of  kinetic  energy, 

(3U/3z)2,  is  dissipated  by  the  buoyancy  and  less  by  the  turbulence 

generated.  An  increasing  Ri,  therefore,  means  a  lower  level  of 

turbulence  and  a  lower  A  or  D  .  The  Ri  in  Equation  23  is 

z  z 

called  a  "local”  Richardson  number  or  gradient  Richardson  number 
as  opposed  to  a  "bulk"  Richardson  number  applied  to  the  whole 
water  column. 


The  argument  that  a  high  Ri  should  suppress  A  and  D  and 

z  z 

no  stratification  should  have  no  effect  along  with  some  dimen¬ 
sional  arguments  led  Munk  and  Anderson  (1948)  to  the  Richardson 
number  function  of 

F( Ri )  =  (1  +  a  Ri)b  (24) 

where  they  deduced  for  vertical  momentum  transport  a=10  and 
b=-l/2  and  for  density  a=10/3  and  b=-3/2.  Since  the  exponent, 
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b,  is  negative  then  in  Equation  22  a  higher  Ri  leads  to  a 
lower  A  .  The  form  of  the  Richardson  number  function  given  in 
Equation  24  appears  in  many  studies.  It  tends  to  have  a 
theoretical  basis  as  shown  by  the  heuristic  derivations  given 
by  Officer  (1976). 

Another  form  of  the  Richardson  number  function  used  by 

Leendertse  and  Liu  (1975)  is 

F( Ri )  =  e-rRl  (25) 

which  is  an  exponential  form  that  has  a  value  for  negative  Ri 

while  Equation  24  can  degenerate  to  zero.  A  negative  Ri  exists 

where  an  inverse  density  gradient  occurs  and  mixing  takes  place 

by  turnover  or  by  vertical  convection  at  a  very  rapid  rate.  It 

can  be  handled  in  numerical  computations  by  using  an  algorithm 

that  sets  a  large  A  where  the  density  gradient  is  negative. 

For  large  Ri,  both  Equations  24  and  25  can  lead  to  small 

A  ,  the  lower  limit  of  which  must  be  the  molecular  diffusivity 
z 

and  viscosity. 

The  most  complete  evaluation  of  A  and  F(Ri)  for  an  estuary 

Z  0 

has  been  given  by  Pritchard  (1960).  Based  on  data  from  the 
James  River  and  using  mixing  length  arguments,  he  deduced  that 

F( Ri )  =  (1  +  0.276  Ri)-2  (26) 

A  =  8.59*10"3IL  fz2(d-z)2/d3l  (27) 

Z  o  L  L  J 

v/here  Ut  is  the  root -mean- square  (RMS)  tidal  velocity  and  d  is  the 

total  depth.  The  Richardson  number  was  defined  as  (g/p  3p/3z)/ 

(0.70U/h)  where  U  is  the  mean  velocity  over  the  water  column.  The 

function  A  is  shaped  to  the  depth  of  the  water  column  to  be  zero  at 
z  o 


the  surface  and  bottom  and  to  have  a  maximum  at  mid  depth  based  on 
mixing  length  arguments.  This  shape  is  modified  by  F(Ri)  which 
because  of  stratification  reduces  the  mid-depth  value.  The  Ri 
is  computed  from  a  mixture  of  the  local  gradient  and  bulk  or 
water  column  velocity  and  depth.  It  is  clearly  intended  to  be 
a  local  or  depth-dependent  evaluation  of  Az  by  choice  of  the  depth 
f unct ion . 

The  depth-dependent  forms  of  A similar  to  Equations  26 
and  27  were  tried  in  numerical  models  by  Bowden  and  Hamilton 
(1975)  and  by  Elliott  (1976).  Both  reported  numerical  instabil¬ 
ities  that  were  apparently  related  to  the  evaluation  of  Az 
using  a  local  Richardson  number.  Bowden  and  Hamilton  (1975) 
resorted  to  using  a  bulk  Richardson  number  but  with  a  vertical 
shape  function.  More  recently,  Bowden  (1977)  has  discussed  the 
difficulties  of  utilizing  the  local  Richardson  number  in  a  nu¬ 
merical  model  and  suggests  that  there  is  a  theoretical  rationale 
for  the  bulk  Ri.  Examination  of  the  numerical  form  of  the 
Elliott  and  the  Bowden  and  Hamilton  models  suggests  that  the 
instabilities  are  a  result  of  computational  procedures  rather 
than  physical  aspects  of  the  mixing  problem. 

As  shown  in  Part  V,  LARM  is  run  with  a  local  Ri  determined 

from  the  local  velocity  and  density  gradient  using  Equation  24. 

When  computations  are  made  for  an  estuary  with  a  constant  A  , 

z  o 

it  produces  A  that  has  the  depth  distribution  as  given  by 
Pritchard  (1960).  This  suggests  that  the  distribution  of  Az  is 
a  result  of  the  interaction  of  the  mean  velocity  and  salinity 
fields  with  the  vertical  dispersion. 
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Wind  is  important  in  vertical  mixing.  Wind  not  only  increases 
the  mean  velocity  gradient  and  shear,  but  must  also  increase  the 
turbulent  fluctuations  about  the  mean.  It  is  the  fluctuating 
component  over  a  given  period  of  averaging  that  determines  verti¬ 
cal  mixing.  Pritchard  (1960)  assumed  that  the  wind-induced 
turbulence  is  proportional  to  the  orbital  wave  velocity  resulting 
from  the  wind,  and  is  also  affected  by  the  vertical  density 
gradient.  The  wind  speed  function  deduced  by  Pritchard  (1960) 
from  the  James  River  data  is 

A_  =  9.57xl0-3  z(h-z)  H  e“2lrz/L  (28) 

Z°  h  T 

where  H  is  the  wind  wave  height,  T  is  the  wave  period,  L  is  the 

wave  length,  and  the  resulting  A  is  additive  to  Equation  27. 

z  0 

In  the  above  relationship,  the  A  is  seen  to  decrease  exponen- 

Z  0 

tially  with  depth  and  also  to  have  a  parabolic  shape  with  depth. 

The  exponential  decay  is  from  the  wave  orbital  velocity  decay 

with  depth  as  derived  from  elementary  wave  theory.  The  parabolic 

shape  terms  come  from  the  mixing  length  theory  used  to  scale  to 

the  total  water  column  depth.  Ford  (1976)  has  examined  the 

wind  effects  on  A  primarily  for  lakes.  The  exponential  decay 

zo 

is  used  but  there  is  no  shape  function  with  depth.  The  wind  wave 
characteristics,  H,  T,  and  L  are  functions  of  wind  speed,  dura¬ 
tion,  and  fetch.  The  functional  relationships  between  these 
variables  are  known  to  a  limited  extent  for  estuaries.  Accurate 
representation  of  vertical  mixing  on  estuaries  will  require 
performing  wind  wave  analyses  at  least  for  significant  wind  events. 


The  scale  of  the  computational  model  grid  relative  to  the 
waterbody  has  a  significant  effect  on  the  magnitude  of  the  dis¬ 
persion  coefficients.  Increasing  layer  thickness  from  0.5  m  to 
1  m  and  to  2  m  reduces  the  resolution  of  the  computations  and 

the  added  coarseness  may  mean  larger  A,,  if  not  more  uncertain 

Zo 

values.  In  principle,  making  the  grid  smaller  approaches  molec¬ 
ular  scales.  However,  the  turbulence  in  the  waterbody  consists 
of  random  motions  and  still  requires  evaluation  of  the  turbulent 
transport  terms  <u'  w’>  and  <w'  s’>  to  close  the  equations  of 
motion  and  transport  for  solution.  A  unique  experiment  in  model 
scales  was  to  compare  LARM  computations  on  a  very  small  grid  of 
Ax=1.52  m  and  H=7.62  cm  to  a  hydraulic  flume  with  density  under¬ 
flow  that  was  run  for  obviously  laminar  conditions  (Edinger  and 
Buchak,  1979b).  The  numerical  model  did  not  produce  accurate  re¬ 
sults  until  the  dispersion  coefficients  were  reduced  to  molecular 
values,  indicating  that  a  strong  interrelationship  exists  be¬ 
tween  grid  size,  dispersion  coefficients,  and  turbulent  charac¬ 
teristics  of  the  waterbody. 

Longitudinal  Dispersion  Formulations 

The  longitudinal  dispersion  of  momentum  and  constituent  results 
from  the  time  averaging  that  produces  the  advection  of  momentum, 
3UU/3x,  in  the  momentum  balance  and  of  constituent,  3US/3x,  in  the 
transport  balance.  The  average  of  the  product  fluctuations  about 
the  mean  values  are  for  momentum: 

< u '  u * >  =  -  Ax  3U/3x  (29) 
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and  for  constituent: 

<u '  s'  >  =  -  Dx  as/  dx  (30) 

which  relates  the  turbulent  fluctuation  to  the  average  velocity 
U  and  constituent  S.  As  discussed  previously,  there  is  a  rela¬ 
tionship  between  the  time  over  which  the  fluctuation  products 
are  averaged  and  the  time  step  of  numerical  integration  that 
produces  U  and  S. 

Most  formulations  of  longitudinal  dispersion  in  estuaries 
have  been  derived  in  relation  to  one-dimensional  sectionally 
homogenous  estuary  models.  They  often  come  from  analysis  of 
steady  channel  flows  and  are  dimensionally  of  the  form: 

Dx  =  a  h*u*  (31) 

where  h*  and  u*  are  a  characteristic  depth  and  velocity. 

Harleman  (1971)  shows  that  when  the  Taylor  formula  is  converted 
from  pipe  flow  to  a  tidal  case,  the  characteristic  velocity 
translates  from  the  boundary  shear  velocity  to  the  maximum  tidal 
velocity.  Examples  of  Equation  31  are  given  in  Fischer  et  al . 
(1979). 

Numerous  evaluations  of  Dx  have  been  made  for  specific 
estuaries  using  a  one -dimensional  model  of  the  salinity  distribu¬ 
tion  (Officer,  1977).  For  the  Potomac  estuary  values  of  10  to 
20  m2/s  were  found  in  the  vicinity  of  Washington,  D.  C.,  increasing 
to  60  m2/s  toward  the  mouth.  Using  a  dissolved  oxygen  model  on 
the  Delaware,  values  ranged  from  120  to  210  m2/s  over  135  km. 

The  Thames  has  exhibited  values  between  50  and  80  m2 / s  at  low 
river  flows  and  up  to  340  m2/s  at  high  river  flows.  The  values, 


therefore,  range  from  5  to  500  m2/s  and  vary  with  river  flow 
and  salinity  stratification.  The  latter  can  be  complicating 
since  a  large  Dx  required  in  a  one-dimensional  model  in  the 
stratification  region  is  simply  forcing  the  model  to  fit  a 
multilayered  circulat ion, and  a  similarly  large  Dx  would  not  be 
required  for  a  laterally  averaged  model  in  the  same  region. 

There  have  been  few  attempts  to  derive  relationships  for 
the  use  of  A  and  D  in  laterally  averaged  hydrodynamics  except 

X  X 

to  use  relationships  of  the  form  in  Equation  31  as  if  they 
applied  to  each  layer.  Experience  with  LARM  has  been  that  two 
order  of  magnitude  changes  in  A^  and  Dx  have  not  changed  results 
significantly.  This  insensitivity  may  be  due  to  upwind  dif¬ 
ferencing.  A  contributing  factor  for  this  lack  of  sensitivity 
is  that  the  contribution  of  vertical  exchange  to  longitudinal 
mixing  which  is  absorbed  in  Dx  in  one-dimensional  models  is 
explicitly  included  in  the  laterally  averaged  models.  Thus, 
one  of  the  spatial  dimensions  over  which  the  turbulent  transport  is 
averaged  in  the  one-dimensional  model  is  relaxed  in  the  two- 
dimensional  model  and  it  becomes  less  important. 

There  still  may  be  a  scale  relationship  to  the  grid  Ax  for 
Ax  and  Dx  as  discussed  previously  for  the  flume  studies  at  laminar 
flow  conditions.  It  suggests  that  in  each  computational  grid  it 
is  necessary  to  satisfy  the  condition  of 

Dmol  <  D  <  UAx  (32) 

x 

which  dimensionally  states  that  constituent  (including  momentum) 
should  not  be  diffused  out  of  a  cell  faster  than  it  is  advected 
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inward.  For  an  estuary  like  the  Potomac  where  computationally 
Ax  might  be  9300  m  (5  nautical  miles  (nm))  and  the  average  net  non- 
tidal  velocity  is  of  the  order  of  0.02  m/s,  the  above  relationship 
gives  186  m2/s.  The  LARM  computations  use  100  m2/s;  Elliott  (1976) 
used  10  to  1  nr/s  for  similar  simulations.  There  is  more  to 
be  learned  about  the  role  of  Ax  and  in  the  laterally  averaged 
models  as  opposed  to  the  one-dimensional  cases. 
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PART  III:  MODEL  SETUP  AND  TESTING 

Setting  up  the  model  for  a  particular  estuary  requires 
determining  the  cross-sectional  geometry,  describing  the  boundary 
conditions  and  initial  conditions,  and  performing  validation 
computations.  After  these  conditions  are  established,  it  is 
necessary  to  design  simulations  that  determine  how  the  model 
performs  over  a  long  period  of  time.  Often  it  is  necessary  to 
redesign  the  simulation  conditions  in  order  to  compare  different 
results  or  to  obtain  different  forms  of  output  for  post 
processing . 

Model  Geometry 

The  Potomac  River  estuary  is  shown  in  Figure  2.  It  is 
98  nautical  miles  (nm)  from  the  mouth  at  Chesapeake  Bay  to  the 
head  of  tide.  It  is  5  to  8  nm  wide  near  the  mouth  and  maintains 
its  breadth  to  Morgantown  near  nm  45.  For  modeling  purposes  the 
estuary  is  conveniently  divided  into  5-nm  segments  numbered  con¬ 
secutively  from  the  head  of  tide  to  the  mouth  as  shown  in 
Figure  2 . 

The  waterbody  geometry  in  laterally  averaged  waterbody 
dynamics  is  represented  by  the  width  of  the  waterbody  at  each 
depth  at  the  center  of  each  computational  cell.  Cross-sectional 
geometry  of  the  Potomac  River  estuary  is  tabulated  for  every 
nautical  mile  in  Cronin  and  Pritchard  (1975).  Estuarine  widths 
for  the  Potomac  used  in  the  model  are  given  in  Table  1.  Table  1 
also  shows  the  computational  scheme  of  the  model  as  the  combination 
of  I's  and  K's  that  make  up  the  computational  grid. 
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The  scheme  is  20  segments  long  of  5  nm  each  and  up  to 
11  layers  deep  each  of  2  metres.  The  shallow  reach  from  model 
segment  15  to  17  is  shown  in  Figure  2  to  be  Kettle  Bottom  Shoals, 
and  segment  12  to  14  is  the  deep  trench  around  Morgantown. 

A  first  impression  is  that  the  rectangular  array  of  computa¬ 
tional  cells  does  not  represent  the  bottom  topography,  and  in 
particular,  the  bottom  slope.  The  bottom  topography  and  slope 
is,  however,  represented  in  detail  by  the  widths.  Connecting 
a  given  bottom  width,  say  400  m,  which  might  represent  the 
channel  prism  line,  shows  that  the  model  has  a  realistic  bottom 
variation  and  uniform  bottom  slope  toward  the  river  end  regard¬ 
less  of  the  computational  grid. 

Boundary  and  Initial  Conditions 

The  three  time-varying  inputs  are  tide  at  mouth,  time-varying 
vertical  salinity  profiles  at  mouth,  and  freshwater  inflow. 
Relative  to  model  geometry,  the  tide  and  salinity  are  both 
specified  in  segment  22  which  is  outside  of  the  Potomac.  The 
geometry  of  segment  22  does  not  affect  the  computations  at  the 
mouth  or  upestuary.  The  first  computed  velocity  profiles  are 
at  the  right-hand  boundary  of  segment  21  and  represent  the  com¬ 
puted  boundary  velocities. 

The  tide  height  variation  at  the  mouth  is  based  on  data 
provided  in  Rives  (1973).  The  present  simulations  use  a  simple 
sinusoidal  tide  with  a  0.20-metre  amplitude  and  12.45-hr  period. 
The  diurnal  inequality  is  not  considered  in  these  simulations, 
nor  are  the  longer  term  mean  water  level  changes  over  a  few  days 
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by  wind.  The  latter  requires  a  measured  tidal  record  at  the 
mouth.  The  data  input  formats  and  time-varying  data  selector 
of  LARM  are  capable  of  handling  these  records,  but  the  records 
were  not  readily  available  for  these  simulations  (Boicourt, 

1980) . 

Time-varying  vertical  salinity  profiles  are  required  at 
the  boundary.  Hourly  salinity  periods  of  September  4  to  7,  1974, 
and  September  9  to  14,  1974,  for  stations  located  at  10  and  19  nm 
(P10  and  P19)  from  the  mouth.  These  stations  are  at  the  down- 
estuary  end  of  model  segments  19  and  17  respectively  (Figure  2) . 

They  are  too  far  upestuary  to  provide  boundary  salinity  data  but 
can  be  used  for  comparisons  to  model  results.  Longitudinal-vertical 
salinity  profiles  extending  to  RM  02  nm  are  given  in  Elliott  (1976) 
approximately  once  a  month  including  August  and  September  1974. 
Examination  of  the  salinity  profiles  from  one  time  to  another 
and  the  hourly  data  in  Elliott  and  Hendrix  (1976)  shows  only 
small  hourly  salinity  changes  at  P10  over  the  period.  For  purposes 
of  the  test  simulation,  the  RM  02  nm  salinity  profile  given  in 
Elliott  (1976)  is  used  at  the  boundary  and  assumed  constant  over 
time.  The  time-varying  data  selector  developed  for  model  appli¬ 
cations  is  quite  capable  of  handling  time-varying  vertical  salin¬ 
ity  profiles  when  that  boundary  data  is  available. 

Initial  conditions  refer  to  the  values  of  all  computed 
variables  specified  at  interior  noints  from  which  the  computa¬ 
tion  is  started.  The  numerical  time  integration  of  the  laterally 
averaged  waterbody  dynamics  written  in  LARM  and  the  time  stepping 
procedure  is  such  that  one  can  begin  with  zero  velocity  components 


and  constant  constituent  concentration,  or  in  the  case  of  an 
estuary  zero  salinity,  and  let  the  computations  interate  to  a 
time-varying  solution  from  the  boundary  data  alone.  This  was 
not  done  for  the  Potomac  because  it  takes  approximately  30  to 
60  simulated  days  to  achieve  an  internal  distribution  for 
salinity.  However,  it  requires  only  two  or  three  tidal  cycles 
to  establish  periodic  water  surface  profiles  and  velocity 
profiles.  A  solution  to  this  problem  is  to  specify  the  initial 
salinity  distribution  either  from  measurements  or  previous 
computations.  The  profile  in  Elliott  (1976)  for  August  21,  1974, 
was  used  to  give  the  initial  longitudinal  and  vertical  salinity 
profiles.  These  are  given  in  Table  2.  The  transient  build  up 
of  salinity  from  boundary  data  is  examined  more  fully  in  Part  V. 

Freshwater  inflow  to  the  Potomac  estuary  is  from  the  river 
at  Great  Falls  and  numerous  tributaries  along  the  100  miles  of 
estuary.  Although  the  model  derivation  and  code  includes  tribu¬ 
tary  inflows,  they  were  not  used  in  the  test  simulation.  The 
daily  mainstream  flow  of  the  Potomac  River  for  the  period  of 
August  and  September  1974  was  taken  from  Elliott  (197b). 

Choice  of  Time  Step 

Choice  of  the  integration  time  step  depends  on  stability 
limits  of  the  numerical  formulation  and  for  tidal  cases  with 
short-term  variations  it  determines  the  resolution  of  the 
boundary  data.  Too  long  a  time  step  means  that  the  time-varying 
boundary  data  may  be  poorly  represented.  The  I. ARM  model  is 
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implicity  longwave  and  the  time  step  is  not  restricted  by  the 


surface  gravity  wave  criterion.  It  is,  however,  restricted  by 
the  Torrence  condition  At<Ax/U. 

It  is  difficult  to  set  At  for  an  estuary,  as  opposed  to  a 
reservoir,  because  U  is  changing  rapidly  with  time  and  can  be¬ 
come  quite  large  within  the  tidal  cycle.  Measurements  on  the 
Potomac  indicated  U'40  cm/s  max,  and  with  Ax=9266.25  m  (5  nm) , 
this  suggests  a  At  limit  of  386  minutes  or  five  hours.  This  is 
much  too  long  to  resolve  the  boundary  tidal  period  as  would  be 
one  hour  or  even  15  minutes.  Because  of  the  efficiency  of  the 
computation  and  the  need  to  have  good  tidal  averages  of  various 
results,  it  was  decided  to  use  a  ten-minute  (600-second)  time 
step  in  the  computations. 

The  numerical  formulation  does  have  an  internal  wave  limita¬ 
tion  as  Ax/At  >  ZgHAp /p  since  the  horizontal  density  gradient 
component  of  the  horizontal  pressure  gradient  is  lagged  by  a 
time  step  in  computing  surface  elevations  and  velocity  components. 
The  Ap/p  for  the  Potomac,  from  S=0  ppt  to  S=15  ppt  is  0.011,  and 
H  20  m  requiring  Ax/At  >  1.466  and  At  <  6317  sec  or  105  minutes 
The  internal  wave  criterion  is  more  restrictive  than  the  Torrence 
cond i t ion . 

Some  runs  were  made  at  a  At  which  divided  the  tidal  period 
of  12.45  hours  into  an  even  number  of  increments.  For  example, 
a  At=448.2  seconds  gives  exactly  100  iterations  per  tidal  cycle 
and  allows  the  tidal  boundary  data  to  be  resolved  at  approxi¬ 
mately  five-minute  intervals. 
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Validation  of  the  computational  code  determines  if  it  gives 
correct  answers  for  known  situations.  It  is  a  step  which  checks 
the  numerical  computations  by  using  relatively  simple  combinations 
of  boundary  conditions.  Validation  cases  are  set  up  to  test  for 
instabilities,  computational  errors,  and  errors  in  accumulation 
or  loss  of  volume  and  salt.  These  would  be  errors  related  to 
formulation  and  coding.  The  estuarine  boundary  condition  of 
LARM  was  run  for  the  cases  of  constant  river  inflow  at  constant 
downestuary  elevation,  constant  salinity  at  the  downestuary 
boundary  with  no  freshwater  inflow,  and  a  sinusoidal  downestuary 
elevat ion . 

The  case  of  constant  river  inflow  and  constant  downestuary 
elevation  results  in  an  upestuary  surface  slope  and  a  computed 
outflow  at  the  mouth.  The  crucial  tests  are  that  the  outflow 
should  eventually  become  equal  to  the  inflow  over  time,  and 
there  should  be  no  volume  loss  of  water  over  time  once  the  sur¬ 
face  slope  is  established.  These  tests  were  all  satisfied. 

The  case  of  constant  salinity  at  the  downestuary  boundary 
tests  the  ability  of  the  computed  boundary  condition  to  convect 
salinity  into  the  estuary  through  the  horizontal  density  gradient. 
It  also  tests  representation  of  the  salt  front  as  it  moves  up¬ 
estuary  due  to  buoyant  convection.  The  model  showed  the  influx 
of  boundary  salts  initially  over  the  total  water  column.  As 
time  progressed,  the  influx  of  boundary  salts  was  from  the 
deeper  portions  of  the  water  column  and  an  outflow  became  estab¬ 
lished  in  the  top  portion  of  column  as  would  be  expected. 
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The  sinusoidal  elevation  at  the  downestuary  boundary  rep¬ 
resents  the  tidal  forcing  function.  In  certain  spatial  dif¬ 
ferencing  forms  of  the  basic  momentum  and  continuity  equations 
it  can  lead  to  producing  higher  frequency  and  shorter  wave¬ 
length  oscillations.  It  can  also  lead  to  alternating  solutions 
on  successive  time  steps.  None  of  these  problems  were  encountered 
using  the  estuarine  boundary  condition  in  the  LARM  code.  The 
tidal  wave  computed  by  the  model  propagated  up  the  estuary  at 
the  surface  gravity  wave  speed.  The  computed  flows  in  and  out 
of  the  mouth  of  the  estuary  became  periodic  within  a  few  tidal 
cycles  and  did  not  show  any  differences  from  one  tidal  cycle 
to  the  next. 

Simulation  Results 

The  first  results  of  interest  from  the  simulations  are  the 
spatial  and  temporal  variation  of  tide  height  and  the  pattern 
of  circulation  throughout  a  tidal  cycle.  The  simulation  was 
run  for  ten  days  to  remove  effects  of  initialization  and  detailed 
hourly  results  were  examined  for  the  last  tidal  cycle.  It  was 
determined  that  initialization  effects  were  overcome  and  station¬ 
ary  conditions  reached  by  comparing  velocity  components  at  a 
given  location  from  one  tidal  cycle  to  the  next. 

The  velocities  were  started  from  zero  and  essentially 
established  a  repetitive  pattern  within  two  tidal  cycles. 

After  that  small  adjustments  of  tenths  of  cm/sec  took  place 
in  response  to  time-varying  freshwater  inflow  and  gravitational 
circulation  related  to  changes  in  salinity. 


Salinity  was  initialized  at  the  values  in  Table  2.  The 
salinity  distributions  at  the  end  of  ten  days  are  shown  in 
Table  3.  Compared  to  the  initial  salinity  distributions,  the 
salinity  after  ten  days  shows  a  very  large  adjustment  from  model 
segment  10  at  river  mile  58  upestuary  to  the  head  of  tide.  This 
computed  change  in  salinity  does  not  appear  to  follow  the  smaller 
salinity  changes  found  in  Elliott  (1976)  between  August  and 
September  1974,  for  the  upestuary  stations. 

The  modes  of  circulation  over  a  tidal  cycle — can  be  examined 
simultaneously  as  the  water  surface  elevation  and  velocity 
vectors.  These  are  shown  for  every  two  hours  through  the  tidal 
cycle  in  Figures  3  to  9.  The  velocity  vectors  are  exagger¬ 
ated  in  the  vertical  by  the  horizontal  to  vertical  scale  of  the 
depth  profile.  An  example  of  the  exaggeration  is  found  near 
110  km.  The  vectors  show  a  vertical  displacement  of  about 
2  metres  over  four  hours  at  the  same  time  that  upestuary  hori¬ 
zontal  motion  is  about  5  km. 

The  water  surface  elevation  for  the  first  hour  of  the  tidal 
cycle  shows  a  high  tide  at  the  mouth  of  the  estuary  simultaneously 
with  a  low  tide  at  the  head  of  the  estuary.  It  is  a  wave  with 
a  node  point  near  RM  50  which  is  at  a  sharp  bend  in  the  river. 

The  circulation  follows  the  water  surface  elevation  changes 
while  it  is  also  responding  to  the  horizontal  density  gradient. 

In  Figures  3  to  9  as  the  circulation  is  followed  from  hour 
to  hour,  it  is  seen  that  the  upper  end  of  the  estuary  continues 
to  flood  while  the  lower  portion  is  ebbing  as  the  circulation 
responds  to  the  surface  wave  (day  10.1181). 
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Figure  3.  Water  Surface  Elevation  and 
Circulation  Pattern  at  0.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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Figure  4.  Water  Surface  Elevation  and 
Circulation  Pattern  at  2.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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Figure  5.  Water  Surface  Elevation  and 
Circulation  Pattern  at  4.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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Figure  6.  Water  Surface  Elevation  and 
Circulation  Pattern  at  6.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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Figure  7.  Water  Surface  Elevation  and 
Circulation  Pattern  at  8.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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Figure  8.  Water  Surface  Elevation  and 
Circulation  Pattern  at  10.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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Figure  9.  Water  Surface  Elevation  and 
Circulation  Pattern  at  12.8  Hours 
after  High  Tide  at  the 
Estuary  Mouth 
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There  is  some  indication  of  two-layered  flow  toward  the  end 
of  the  tidal  cycle  (day  time  10.3681)  which  shows  that  the 
flooding  tide  at  the  estuary  mouth  is  mostly  in  the  bottom 
layers . 

There  is  a  tendency  for  the  magnitudes  of  the  velocity 
vectors  at  a  given  time  in  the  tidal  cycle  to  follow  the  geom¬ 
etry  of  the  estuary.  The  strongest  flows  are  at  the  constricted 
cross  section  near  model  segment  11  which  is  also  one  of  the 
more  shallow  portions  of  the  estuary.  Velocities  are  smaller 
throughout  the  tidal  cycle  in  the  lower  portion  of  the  estuary 
because  of  its  greater  cross-sectional  geometry  (Table  1). 

Overall  flushing  and  transport  through  the  estuary  is 
eventually  related  to  the  tidally  averaged  circulation.  The 
computed  tidally  averaged  circulation  for  the  Potomac  is  shown 
in  Figure  10.  Here  the  vector  lengths  represent  displacement 
over  24  hours.  It  is  a  result  of  averaging  the  velocity  vectors 
at  each  point  over  a  tidal  cycle  at  ten-minute  intervals. 

The  computed  boundary  velocities  are  at  the  mouth  of  the 
estuary  near  190  km.  From  Figure  2,  this  is  at  the  boundary 
where  the  Potomac  enters  Chesapeake  Bay.  The  boundary  veloc¬ 
ities  are  computed  from  the  specified  tide  height  in  the  bay 
and  the  simultaneously  computed  tide  height  for  the  first  model 
segment  in  from  the  bay  (segment  21).  The  velocities  at  the 
mouth  show  a  tidally  averaged  profile  of  intrusion  in  the 
bottom  layers  and  an  outflow  in  the  surface  layers.  This 
pattern  follows  upestuary  where  at  each  cross  section  there  is 
a  distinct  depth  at  which  the  vectors  change  from  having  an 
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upestuary  component  to  a  dcwnestuary  component.  Overall,  the 
estuary  tends  to  circulate  in  two  parts:  an  upestuary  part 
above  80  km  at  the  narrows  and  a  downestuary  part.  There  is  a 
net  flow  between  the  two  parts  represented  by  the  bottom  water 
intrusion  near  80  km. 

The  tidally  averaged  circulation  shown  in  Figure  10 
tends  to  have  a  circular  pattern  near  the  estuarine  mouth. 

There  is  a  slight  downflow  from  the  surface  layer  to  the  bottom 
layers  in  the  last  segment  before  the  boundary  outflows.  Since 
the  computation  of  the  outflow  boundary  velocities  does  not 
include  the  horizontal  advection  or  dispersion  of  momentum, 
the  circulation  may  result  from  dissipation  of  momentum  from 
the  interior  flow  at  the  boundary.  The  magnitude  of  the  rotary 
circulation  is  not  as  large  as  indicated  in  Figure  10  because 
of  the  exaggerated  vertical  velocity. 
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PART  IV:  PREDICTED  AND  OBSERVED  TIDAL  DYNAMICS 


The  availability  of  survey  data  on  the  Potomac  River  estuary 
allows  quantitative  comparisons  of  predicted  model  results  and 
observed  tide  height  and  velocity  conditions.  Tide  height  in¬ 
formation  is  available  in  some  form  along  most  estuaries.  It 
may  be  time  and  height  of  low  and  high  water  forecasts  as  a  min¬ 
imum  or  detailed  tide  height  measurements  over  a  specific  period 
of  time.  The  difference  is  that  the  latter  includes  wind  and 
freshwater  runoff  effects  while  the  former  includes  only  soli- 
lunar  effects.  Wind  and  freshwater  runoff  influences  the  mean 
water  level  over  a  number  of  days,  and  their  effect  on  the 
downestuary  boundary  are  important  for  proper  computation  of 
boundary  inflows  and  outflows. 

Detailed  velocity  profiles  are  not  routinely  obtained  for 
estuaries.  Special  studies  of  the  Potomac  were  conducted  by 
the  Chesapeake  Bay  Institute  and  reported  in  Elliott  and 
Hendrix  (1976)  for  two  periods,  September  4  to  September  7,  1974, 
and  September  9  to  September  14,  1974,  at  two  stations.  One 
station  was  10  nm  from  the  mouth  and  the  other  was  19  nm  from 
the  mouth.  Velocities  were  measured  at  1-metre  intervals  hourly. 
These  velocity  profiles  are  used  for  detailed  comparison  to  the 
model  results. 

Vertical  dispersion  coefficients  and  vertical  velocity 
components  were  not  computed  for  the  Potomac  data.  They  are, 
however,  an  important  part  of  the  model  since  they  are  a  major 
transport  mechanism.  The  model  results  are  compared 


qualitatively  with  the  properties  of  vertical  dispersion  coeffi¬ 
cients  and  vertical  velocity  components  inferred  by  Pritchard 
(1967)  from  the  James  River  studies.  The  James  River  studies 
are  also  summarized  in  Officer  (1976). 

Tide  Heights 

The  predicted  model  tide  heights  over  time  for  various 
stations  along  the  Potomac  are  given  in  Figure  11.  The  tide 
is  driven  at  the  mouth  (186.25  km)  with  an  amplitude  of  0.2  m. 
Upestuary  at  28.73  km,  the  tide  range  tends  to  increase  and 
there  is  a  phase  shift  of  high  tide  to  occur  later  than  at  the 
mouth.  The  range  and  phase  shift  of  the  observed  tides  have  been 
summarized  for  the  Potomac  by  Rives  (1973),  and  can  be  used  for 
comparison  to  the  computations. 

Tidal  ranges  are  compared  in  Figure  12.  The  model  results 
are  given  for  various  values  of  bottom  friction  since  this  is 
the  major  variable  that  influences  tidal  amplitude  in  the  model. 

The  observed  tide  range  shows  a  local  maximum  upestuary  from  the 
mouth,  a  minimum  near  model  segment  11  (Figure  2),  and  increase 
to  the  head  of  tide.  The  observed  tide  range  minimum  coincides 
with  the  nodal  point  of  the  standing  tide  wave. 

The  computed  results  show  that  the  model  runs  satisfactorily 

over  a  wide  range  of  friction  factors  with  no  numerical  difficulties 

Numerical  schemes  that  lag  many  terms  often  fail  at  low  friction 

(high  Chezy  C)  since  this  term  often  maintains  numerical 

stability.  The  change  in  tide  range  from  the  mouth  to  the  head 

of  tide  is  reproduced  by  the  model  for  a  C=60  m^/s.  The  minimum 

x 

at  90  km  is  reproduced  for  a  C=40  m“/s. 
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Figure  11.  Time-varying  Tide  at  Selected  Stations 
along  the  Potomac  River  Estuary  Computed 
for  a  Sinusoidal  Tide  at  the  Mouth 
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The  model  could  be  tuned  to  tide  range  and  phase  by  varying 
the  friction  coefficient  longitudinally  along  the  model. 

Rives  (1973),  however,  shows  with  a  one- dimensional  tidal  model 
with  detailed  spatial  averaging  that  the  Potomac  estuary  tidal 
range  minimum  is  reproduced  for  a  spatially  uniform  friction 
coefficient.  Varying  the  friction  coefficient  spatially  in  a 
coarse  spaced  grid  as  used  in  the  present  simulation  thus  appears 
partially  to  be  a  surrogate  for  spatial  detail.  There  are,  how¬ 
ever  ,  changes  in  cross-sectional  geometry  for  which  changes  in 
the  friction  coefficient  is  justified. 

The  phase  lag  of  the  tide  height  from  the  mouth  of  the 
estuary  can  be  compared  for  observed  and  computed  conditions. 

It  is  given  for  various  friction  coefficients  in  Figure  13. 

There  is  a  tendency  for  the  model  to  underestimate  the  phase 
lag  at  the  head  of  the  estuary.  There  is  little  variation  in 
the  computed  phase  lag  with  friction  coefficient  since  phase  lag 
mostly  depends  on  gravity  wave  propagation.  Propagation  of  the 
tide  wave  up  the  estaury  is  primarily  governed  by  the  gravity 
wave  speed  (/gH)  which  is  solely  function  of  depth.  The  computed 
wave  would  be  based  on  the  geometry  given  in  Table  1  and  this 
may  vary  from  the  real  geometry  related  to  wave  propagation.  The 
difference  in  observed  and  computed  tide  phase,  based  on  the 
results  of  Rives  (1973),  may  be  due  to  the  geometry  resolution 
in  the  model . 

Velocity  Distributions 

The  real  test  of  the  model  is  the  degree  to  which  it  re¬ 


produces  observed  temporal  variations  in  velocity.  A  valid 


Phase  Lag,  Hours  from  Mouth 


Figure  13.  Comparison  of  Observed  and  Predicted  Tidal 
Phase  along  the  Potomac  Estuary 
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comparison  depends  on  how  well  the  boundary  conditions  during 
the  period  of  observation  are  represented  in  the  model.  Conditions 
during  the  Chesapeake  Bay  Institute  (CBI)  September  1974  surveys  are 
given  in  Figure  14  including  the  residual  or  long-term  change  in 
tide  height  at  the  estuary  mouth,  and  the  wind  shear  stress  along 
and  across  the  estuary. 

The  detailed  tide  and  wind  record  for  the  period  of  observa¬ 
tions  used  in  Figure  15  are  no  longer  readily  available 
(Boicourt,  1980).  Lacking  the  actual  data,  the  model  is  run  for 
the  simpler  condition  of  a  sinusoidal  tide  with  no  mean  elevation 
change  and  no  wind.  The  latter  conditions  appear  to  be  satisfied 
for  the  first  18  hours  of  the  September  9  to  14,  1974,  survey 
period  shown  in  Figure  14.  The  hourly  velocity  and  salinity 
profiles  at  station  P19  are  used  for  showing  the  comparison  to 
the  model  results.  The  results  at  station  P10  were  similar  to 
P19.  Deviations  from  the  model  conditions  are  examined  as  part 
of  a  sensitivity  analysis. 

The  computed  and  observed  velocities  at  2-metre  depths  are 
compared  for  survey  station  P19  in  Figure  15.  Since  a  sinusoidal 
tidal  height  at  the  mouth  is  used  in  the  model,  the  relationship 
between  the  tide  phase  and  the  day-hour  is  not  known  precisely. 

The  computed  velocities  at  one  depth  were  shifted  in  phase  until 
they  coincided  with  the  observed  velocities.  The  sinusoidal  tidal 
boundary  in  the  model  does  not  include  any  diurnal  inequality 
which  might  have  been  in  the  real  tidal  record.  The  computed 
velocities  reproduce  the  maximum  and  minimum  velocities  at  all 


Figure  14.  Tide  and  Wind  Conditions  during  Study  Period 
E 7  (cm),  the  Residual  Elevation  at  Colonial 
“  Beach;  Wj  (m*S'l)  Downestuary  Wind; 

W~  (m*S*l)  Cross-Estuary  Wind,  from  NE  to  SW. 
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depths  and  except  for  the  surface  velocities,  the  model  repro¬ 
duces  the  shape  of  the  observations.  The  surface  velocities  are 
distorted  by  wind. 

Important  characteristics  of  estuarine  tidal  velocities  are 
shown  in  the  data  and  reproduced  by  the  model.  In  a  stratified 
estuary  exhibiting  a  two -layered  circulation  like  the  Potomac, 
the  tendency  of  the  velocity  is  to  have  a  stronger  ebb  than 
flood  velocity  in  the  surface  layers  and  to  have  a  stronger 
flood  than  ebb  velocity  in  the  bottom  layers.  This  character¬ 
istic  is  shown  in  the  data  presented  in  Figure  15  and  is  repro¬ 
duced  by  the  model. 

Detailed  comparisons  of  the  vertical  velocity  profiles  are 
given  in  Figure  16.  Model  profiles  were  available  every 
20  minutes  and  observed  profiles  available  every  hour.  The  phase 
shift  of  the  computed  profiles  of  Figure  16  was  nine  minutes, 
thus  the  observed  and  computed  profiles  do  not  coincide  exactly. 
The  model  profiles  produce  the  magnitude  of  the  observed  veloc¬ 
ities  throughout  the  tidal  cycle  and  reproduce  important  profile 
features.  Beginning  toward  mid-tide  at  0800  hours,  the  bottom 
velocities  begin  to  flood  while  the  surface  velocities  are  still 
ebbing.  This  continues  until  the  full  flood  tide  velocity  pro¬ 
file  is  developed  and  then  the  bottom  velocities  begin  to  ebb 
before  the  surface  velocities. 

The  tidally  averaged  velocities  are  a  fraction  of  the  ebb 
and  flood  maximum  velocity,  but  the  net  circulation  and  flushing 
of  the  estuary  is  determined  by  the  tidally  averaged  velocity. 
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Figure  16.  Computed  and  Observed  Velocity  Profiles 
at  Station  P19  Potomac  Estuary 
for  September  9-10,  1974,  Taken 
Hourly  through  One  Tidal  Cycle 
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The  tidally  averaged  and  ebb  and  flood  maximum  velocities  are 
shown  in  Figure  17  along  with  the  tidally  averaged  salinity 
profile.  The  tidally  averaged  velocity  profile  is  the  classical 
case  of  a  surface  layer  outflow  and  a  bottom  layer  inflow.  The 
model  reproduces  the  depth  of  the  velocity  interface  and  the 
bottom  layer  velocities.  The  surface  layer  velocities  deviate 
from  each  other  because  of  wind  effects. 

The  observed  and  computed  maximum  ebb  and  flood  velocities 
agree  except  in  the  lower  portion  of  the  water  column.  The  ob¬ 
servations  show  a  flood  tide  maximum  with  depth  in  the  bottom 
layer  that  has  also  been  observed  by  Pritchard  (1967)  which  is 
not  reproduced  by  the  model. 

Salinity  Distributions 

The  tidally  averaged  salinity  profile  is  shown  in  Tables  2 
and  3.  The  salinity  profile  changes  very  little  throughout  the  tidal 
cycle  at  stations  close  to  the  estuary  mouth  unless  there  are 
changes  in  the  boundary  salinity.  The  computed  salinities 
would  be  expected  to  match  the  observed  salinities  because 
observed  salinities  at  an  earlier  date  were  used  to  initialize 
the  model  computations.  Unfortunately,  there  are  no  detailed 
salinity  profiles  further  upestuary  during  the  CBI  surveys. 

It  is  expected  that  salinities  in  the  upper  reaches  of  the 
estuary  are  more  variable  than  in  the  lower  estuary.  In  the 
upper  estuary,  the  salinity  is  an  order  of  magnitude  lower  and 
a  small  variation  is  a  greater  percent  change.  In  the  upper 
estuary,  salinity  is  controlled  more  by  freshwater  input  than 
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Figure  17.  Computed  and  Observed  Tidal  Averaged 
Velocity  Profile  and  Maximum  Ebb  and 
Flood  Velocities  at  Station  P19 
Potomac  Estuary  for  September  9-10,  1974 
Survey  Period 
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tide  and  salinity  at  the  mouth.  Except  for  comparing  changes 
in  salinity  profiles  over  a  long  period  of  time  as  in  Tables  2 
and  3  from  initial  conditions  to  stationary  conditions,  there 
was  no  detailed  profile  station  of  salinity  and  current  in  the 
upper  estuary  for  more  detailed  comparison  of  salinity  variations. 

Vertical  Velocity  Components 

The  computed  ebb,  flood,  and  tidally  averaged  vertical 
velocity  component  at  RM  10  is  given  in  Figure  18.  The  tidal 
average  is  an  order  of  magnitude  smaller  than  the  maximum  ebb 
or  flood  values,  and  its  consistency  in  comparison  to  observations 
is  a  measure  of  its  accuracy.  The  tidally  averaged  vertical 
velocity  component  is  about  1/200  of  the  horizontal  velocity 
component.  The  vertical  transport  of  mass  and  constituent  re¬ 
lated  to  this  low  velocity  is,  however,  large  in  the  finite  dif¬ 
ference  model  because  of  the  larger  area  through  which  it  is 
transporting.  In  a  model  segment,  the  ratio  of  bottom  area  to 
end  area  is  as  Ax/H  and  for  the  Potomac  simulation  this  ratio  is 
2500,  thus  the  vertical  transport  can  be  as  large  as  the  hori¬ 
zontal  transport. 

The  vertical  velocity  component  cannot  be  measured  directly. 
Pritchard  (1967)  has,  however,  computed  a  vertical  velocity  com¬ 
ponent  profile  in  the  James  River  using  up-  and  downestuary 
velocity  transects  and  salinity  profiles  over  time  and  evaluating 
numerically  various  terms  in  the  momentum  and  salt  balance. 

These  calculations  showed  that  for  a  partially  mixed  estuary 
the  tidally  averaged  vertical  velocity  component  should  be  upward 


Figure  18.  Computed  Vertical  Velocity  Component 
Profiles  as  a  Tidal  Average  and 
Maximum  and  Minimum  for 
Station  P10 
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from  the  bottom  to  the  surface  layers,  should  be  zero  at  the 
bottom,  and  the  surface  should  have  a  maximum  in  the  bottom 
layers.  Figure  18  shows  that  the  vertical  velocity  component 
profile  produced  by  the  model  has  the  expected  properties.  The 
magnitude  of  the  vertical  velocity  component  is  similar  to  that 
found  by  Pritchard  (1967). 

Vertical  Dispersion  Coefficients 

The  theoretical  aspects  of  the  vertical  dispersion  relation¬ 
ship  in  numerical  models  for  estuaries  and  other  waterbodies  was 
discussed  in  Part  II.  The  Munk  and  Anderson  (1948)  Richardson 
number  dependency  of  the  vertical  dispersion  coefficient  on 
stratification  is  used  for  the  present  simulations.  In  order 
to  examine  the  dependence  of  vertical  dispersion  on  the  Richardson 
number  alone,  A  in  Equation  22  was  taken  as  a  constant  of 

Z  0 

_ 

6*10  m2/s  for  both  momentum  and  salt. 

Increasing  the  Richardson  number  in  Equation  24  cither 

by  increasing  stratification  or  decreasing  velocity  shear, 

decreases  A  .  A  lower  value  is  set  for  A  as  the  molecular  value 
z  z 

of  momentum  and  salt  diffusion.  Decreasing  Ri  increases  A  to¬ 
ward  an  upper  limit.  When  Ri  is  <0.10  or  5,/<)z  is  negative  due 
to  vertical  density  instability,  A  is  allowed  t  ,  approach  a 
maximum  value  of  H 2  /  ( 2  r,  t  )  where  H  is  the  layer  thickness  for  'he 
finite  difference  computation.  This  upper  limit  asset  s  tu.it 
there  is  no  numerical  instability  due  to  the  Time  .  apg<.  i  v  i  'i! 
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gradient  terms. 


The  Richardson  number  enters  the  dynamic  computations  as  a 


local  value.  It  is  computed  at  each  layer  for  each  time  step 
and  is  used  to  compute  a  local  Az-  The  vertical  dispersion 
coefficient  is  thus  being  computed  continuously  with  the  local 
velocity  and  density  gradient  on  which  it  depends.  The  varia¬ 
tion  of  the  vertical  dispersion  coefficient  over  depth  and 
throughout  a  tidal  cycle  is  given  in  Figure  19.  It  is  computed 
using  the  velocity  profiles  shown  for  the  same  times  in  Figure  16 
and  for  the  salinity  profile  given  in  Figure  17.  The  salinity 
profile  at  staion  P19  does  not  change  substantially  with  time. 

Most  of  the  change  in  the  kr  profile  throughout  the  tidal  cycle 

z 

is  due  to  changes  in  the  vertical  velocity  gradient  over  time. 

Pritchard  (1967)  gives  an  estimate  of  the  vertical  distribu¬ 
tion  of  vertical  dispersion  as  derived  from  the  James  River  data. 


It  shows  A_  being  minimum  or  near  zero  at  the  surface  and  bottom, 
and  a  minimum  at  mid-depth  near  the  region  of  maximum  salinity 
stratification.  The  A  computed  from  data  by  Pritchard  (1967) 

Sj 

had  maxima  in  the  mixed  surface  layers  and  bottom  layers.  The 
A?  computed  by  the  model  shown  in  Figure  1  exhibits  all  of 
the  observed  properties  throughout  the  tidal  cycle.  The  maximum 
value  of  6*10  m‘/s  is  similar  to  that  computed  from  data  by 
Pritchard  <  1967V  During  most  of  the  tidal  cycle  the  model  A 

z 

shows  neper  and  lower  layer  maxima  as  found  from  data. 


Tii'-  A_  <1 1  it  r  i  but  ion  at  hour  7  and  6  into  the  tidal  e  'to 
ihov.-s  n.i  i;,  i  d  --  e  ■  nth  minimum.  The  velocity  pro  f  ’  '  es  in 

;  •  i  -  t  ■  *  i  1  '■  ■  1  ;  ‘  *  '  I  x  •  V>1  r  •  .!  1  a”''  ’■  ' 
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is  beginning  to  reverse  in  direction  giving  maximum  shear  and 

a  low  Richardson  number.  The  result  is  an  uniform  Az  near  the 

A  value  implying  a  potential  for  maximum  mixing  as  the  tide 
z  o 

is  turning  in  the  bottom  layers. 

The  minimum  in  A  throughout  the  tidal  cycle  is  at  the  depth 
z 

of  maximum  salinity  gradient  which  does  not  substantially  vary 

over  time.  The  magnitude  and  shape  of  the  minimum  is,  therefore, 

determined  primarily  by  a  small  velocity  gradient  at  mid-depth. 

The  individual  velocity  profiles  given  in  Figure  16  shows  that 

there  is  little  or  no  velocity  gradient  through  the  region  of 

the  A  minimum, 
z 

The  Az  Richardson  number  dependency  coupled  with  the  numer¬ 
ical  hydrodynamics  in  a  straightforward  manner  produces  Az 
variations  that  have  been  found  in  the  field.  It  has  not  been 

necessary  to  introduce  spatial  functions  on  A  as  proposed  by 

z 

Blumberg  (1975)  to  obtain  known  depth  variations  in  The 

spatial  functions  appear  to  be  required  for  limited  formulations 
of  the  basic  hydrodynamics. 
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PART  V:  SENSITIVITY  ANALYSES 


There  are  numerous  properties  of  the  model  and  simulation 
which  are  not  illustrated  by  the  comparisons  to  data.  These 
are  the  time  required  to  develop  the  salinity  field  from  zero 
salinity  conditions,  the  effects  of  changes  in  the  vertical 
mixing  coefficients,  and  the  effects  of  wind  on  velocity 
simulat ions . 

Initializing  Salinity 

One  way  to  initialize  the  salinity  field  is  to  begin  with 
a  boundary  salinity  throughout  the  waterbody  and  let  the  salinity 
field  and  flow  field  build  up  together  over  time.  The  profiles 
of  salinity  at  different  times  of  such  a  simulation  are  given 
in  Table  4.  These  simulations  had  the  boundary  salinity 
profile  shown  for  model  segment  22  and  sinusoidal  tide.  It  shows 
the  development  of  salinity  stratification  up  the  estuary  and 
the  gradual  movement  of  the  salt  front  with  time. 

The  time  history  of  surface,  mid-depth,  and  bottom  salin¬ 
ities  over  time  at  model  segment  15  about  halfway  up  the 
estuary  is  given  in  Figure  20.  It  shows  that  after  20  days, 
the  salinity  is  continuing  to  increase  with  time.  Even  after 
20  days  it  is  difficult  to  determine  how  much  longer  is  required 
for  there  to  be  only  tolerable  change  due  to  initialization. 

The  time  required  to  initialize  a  salinity  profile  from 
zero  initial  conditions  is  thought  to  be  the  average  residence 
time  of  the  estuary  based  on  the  net.  non-tidal  flow.  For  the 


gure  20.  Time  History  of  Salinity  Build 
at  40  nm  from  Mouth  of  Estuary  for 
Initiation  from  Zero  Salinity.  For 
Surface  Layer(S),  Mid-depth(M) ,  and 
Bottom  Layer(B) 


Potomac,  the  volume  of  the  estuary  is  7*109  m3,  and  the  net  non- 

tidal  flow  based  on  a  river  flow  of  50  m3/s  is  2.5*107  m3/day. 

This  suggests  a  residence  time,  hence  an  initialization  time, 

of  280  days.  Initialization  of  large  waterbodies  from  a  zero 

initial  salinity  condition  is  not  computationally  very  efficient. 

Sensitivity  to  A 
_ _ _ z 

The  vertical  dispersion  coefficients  and  velocity  and  salin¬ 
ity  fields  are  related  through  Equation  24.  It  is  of  interest 
to  determine  the  influence  of  the  magnitude  of  A  on  the  velocity 

Z  o 

and  salinity  profiles.  Simulations  were  made  with  a  large  of 

z  o 

10  cm2/s,  and  a  low  A  of  0.6  cm2/s.  The  model  is  nominally 

z  o 

run  with  A  of  6  cm2/s.  The  simulations  allow  comparison  of 
z  o 

conditions  over  a  factor  of  15.  The  resulting  tidal  averaged 

velocity  profiles  for  the  high,  low,  and  nominal  A  are  given 

in  Figure  21  for  river  mile  P19  and  P10.  It  is  seen  that  a 

factor  of  15  in  A  has  little  effect  on  the  overall  profile. 

zo  * 

Increasing  A  deepens  the  upper  layer  and  reduces  the  surface 
z  o 

velocity.  For  the  velocity  profile,  a  greater  A  implies 

z0 

greater  vertical  transport  of  momentum,  hence  a  thickening  and 

smoothing  of  the  profiles.  The  changes  in  the  salinity  profile 

were  insignificant.  The  relationships  hold  at  P10  close  to  the 

ocean  boundary  as  well  as  the  P19  further  upestuary. 

The  opposite  effect  is  found  by  Wang  and  Kravitz  (1980). 

Using  a  constant  A  with  no  Richardson  number  condition  and 

z 

changing  it  from  0.1  cm2/s  to  0.5  cm2/s,  they  computed  that  the 
upper  layer  thickened  slightly  but  there  was  an  increase  in 
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surface  and  lower  layer  velocities.  They  attributed  this  to  a 
greater  intrusion  flow  because  of  the  more  uniform  salinity  and 
density  gradient  resulting  from  greater  mixing.  The  results 
may  be  erroneous  because  of  ignoring  the  effects  of  stratifi¬ 
cation  on  A  ,  hence  overestimating  it  relative  to  the  remainder 
of  the  profile,  and  because  of  increased  sensitivity  to  changes 
in  low  values  of  A  in  comparison  to  the  Pritchard  (1967) 

Z  0 

estimate  of  6  cm2/s. 

Effects  of  Wind 

One  important  driving  force  for  most  estuaries  is  wind. 

The  horizontal  density  gradients  due  to  ocean  salt  and  fresh¬ 
water  inflow  and  the  solilunar  tide  at  the  mouth  establish  the 
long-term  circulation  with  a  period  of  two  to  four  weeks.  These 
forces  produce  the  classical  net  circulation  of  bottom  intrusion 
of  salt  water  which  mixes  upward  to  a  less  saline  surface  outflow. 
On  a  more  short-term  basis  Elliott  (1976)  has  shown  from  the 
analysis  of  current  meter  data  that  the  classical  circulation 
pattern  is  disturbed  by  wind  with  a  period  of  four  or  five  days 
for  the  Potomac. 

Two  forces  in  the  model  are  influenced  by  wind.  One  is  the 
surface  wind  stress  along  the  axis  of  the  model.  The  second  is 
the  mean  elevation  at  the  ocean  tidal  boundary  which  must  be 
determined  from  measured  tidal  elevations.  If  the  surface  wind 
stress  is  used  independent  of  the  boundary  elevation  changes 
and  with  a  fixed  mean  elevation,  the  free  water  surface  is  hinged 
at  this  point.  The  inflow  and  outflow  across  the  ocean  boundary 
due  to  surface  slope  is  then  dictated  by  wind  conditions  in  the 
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estuary  alone.  A  fixed  end  condition  can  result  in  an  artifical 
circulation  within  the  estuary  that  may  appear  computationally 
correct  but  may  not  match  the  real  conditions  over  the  period 
of  study. 

The  mean  tidal  elevation  in  the  mouth  of  the  Potomac  is 
determined  by  the  wind-driven  surface  slope  of  Chesapeake  Bay, 
which  in  turn  is  often  influenced  by  wind-driven  sea  level 
changes  off  Cape  Henry.  The  Potomac  is  subjected  to  the  four 
combinations  of  rising  or  falling  mean  elevation  at  the  mouth 
with  up-  or  downestuary  winds.  These  can  produce  four  differ¬ 
ent  circulation  patterns  in  the  lower  reaches  of  the  estuary. 

It  is  necessary  to  have  measured  tide  heights  and  concurrent 
winds  as  input  data  to  the  model  to  prevent  the  computation 
of  artifical  circulation  conditions. 

Within  the  estuary,  the  tendency  of  a  surface  wind  is  to 
produce  a  downwind  surface  current  with  a  return  bottom  current. 

The  return  bottom  current  is  a  result  of  the  surface  elevation 
set  up  adding  to  the  horizontal  pressure  gradient  in  the  bottom 
layers.  It  is  known  as  a  barotropic  response.  In  an  estuary, 
the  wind  circulation  is  combined  with  the  density  circulation, 

A  downestuary  wind  thus  adds  to  the  downestuary  surface  flow 
and  upestuary  bottom  flow.  An  upestuary  wind  can  produce  com¬ 
plex  circulation  patterns  as  a  result  of  the  surface  wind 
current  opposing  the  surface  outflow  and  the  bottom  wind  current 
opposing  the  bottom  inflow.  By  slowing  down  the  bottom  inflow, 
an  upestuary  wind  can  result  in  a  three-layer  flow,  upestuary  in  the 
top  mid-layer  and  downestuary  in  the  bottom  layer.  These  patterns 
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are  changing  from  day  to  day  due  to  the  strength  and  direction  of 
the  wind  and  changes  in  boundary  elevation  with  time,  and  lead 
to  the  panorama  of  circulation  patterns  found  by  Elliott  (1976) 
from  current  meter  data. 

It  is  necessary  to  perform  wind  simulations  with  the  model 
to  show  that  it  gives  reasonable  results.  The  Potomac  base 
case  simulations  were  made  with  no  wind  and  compared  to  a  sur¬ 
vey  period  of  current  meter  and  salinity  data  that  satisfied  this 
condition.  These  comparisons  were  made  because  the  detailed  sur¬ 
vey  period  hourly  wind  speed  and  tide  height  data  were  not 
available.  One  wind  test  would  have  been  to  run  this  period 
with  and  without  wind.  Another  would  have  been  to  run  the  survey 
period  and  determine  wind  effects  by  appropriate  numerical 
filtering  similar  to  that  used  for  the  current  meters  by  Pritchard 
and  Rives  (1979).  Since  wind  data  are  unavailable  for  either  of 
the  above  simulations,  the  tests  are  run  for  the  artificial  condi¬ 
tions  of:  (1)  no  tide  with  steady  wind;  (2)  with  tide  with  and 
without  wind.  The  first  case  is  run  to  eliminate  the  masking  of 
water  surface  elevation  changes  due  to  wind  by  the  tidal  changes 
in  order  to  examine  the  time  history  of  the  water  surface  set  up. 
The  second  case  is  run  to  determine  the  contribution  of  wind  to  the 
circulation  measured  as  the  difference  in  velocity  with  and  without 
the  wind.  Wind  stress  was  parameterized  in  accordance  with  the 
well  known  quadratic  shear  stress  low. 

Running  the  model  for  a  steady  wind  with  no  tide  should  re¬ 
sult  in  a  water  surface  set  up  in  the  direction  of  the  wind 
with  the  establishment  of  the  surface  current.  Countering  the 


surface  current  is  a  return  barotropic  bottom  current  in  response 
to  the  initial  water  surface  set  up.  The  initial  wind  depresses 
the  water  surface  below  the  final  equilibrium  water  surface 
from  which  it  rebounds  due  to  the  bottom  return  flow.  It  is  a 
classical  waterbody  response  to  wind  that  has  been  used  by 
Wang  and  Kravitz  (1980)  to  test  numerical  equations  of  motion. 

The  water  surface  elevation  over  time  for  locations  at  10, 

45,  and  85  nm  from  the  mouth  are  shown  in  Figure  22  for  a  wind 
speed  of  6  m/s.  The  10  nm  location  is  chosen  to  be  near  the 
fixed  elevation  boundary,  the  45  nm  location  is  in  the  narrow 
portion  of  the  estuary,  and  the  75  nm  location  approaches  the 
head  of  the  estuary.  The  water  surface  response  follows  the 
classical  pattern  of  a  single  depression  and  recovery  followed 
by  long  period  oscillations.  The  greater  dip  at  the  head  of  the 
estuary  results  from  the  fixed  elevation  boundary  condition  at 
the  mouth.  The  fixed  elevation  boundary  condition  does  not  allow 
a  barotropic  flow  to  become  immediately  established  and  replace 
the  water  blown  out  of  the  upper  end  of  the  estuary.  It  illus¬ 
trates  the  anomalous  results  that  can  be  obtained  by  mismatching 
the  boundary  elevation  and  surface  wind  condition. 

Simulation  of  the  estuary  with  tide  and  wind  is  more  realistic. 
It  does,  however,  present  the  problem  of  separating  wind-induced 
currents  from  the  larger  tidal  currents.  It  is  done  here  by 
comparing  tidal  averages.  The  base  case  simulations  have  produced 
the  tidal  average  velocity  profile  with  no  wind  (Figure  16). 
Introducing  a  wind  after  the  base  case  has  been  reached  allows  the 
examination  of  the  establishment  of  the  wind-induced  profile. 
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Figure  22.  Wind  Surface  Profiles  Over  Time 
after  Initiation  of  6  m/s  Downestuary 
Wind,  with  No  Tide 
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Figure  23  gives  the  difference  in  velocity  with  and  without 
a  downestuary  wind  of  1.5  m/s.  The  profiles  are  for  successive 
tidal  cycles  for  the  survey  station  of  10  nm  from  the  estuary 
mouth  (Figure  2).  Within  the  first  cycle,  a  barotropic  response 
with  upestuary  flow  is  established  for  the  whole  water  column 
before  the  downestuary  surface  flow  begins.  The  classical  shallow 
surface  flow  in  the  direction  of  the  wind  with  a  deep  bottom 
return  is  established  within  three  tidal  cycles.  There  is  a 
tendency  for  variations  in  the  magnitude  of  the  barotropic  flows 
as  in  the  fourth  and  seventh  tidal  cycles.  The  initial  barotropic 
flow  in  the  first  tidal  cycle  is  related  to  the  initial  water 
surface  depression  shown  in  Figure  22,  while  the  periodic 
changes  in  magnitude  of  the  barotropic  flow  are  related  to  the 
longer  term  dampened  oscillations  also  shown  in  Figure  22. 

The  two  wind  cases  illustrate  the  basic  response  of  the 
model  to  wind  and  show  that  it  produces  expected  results  for 
known  conditions.  Real  winds  vary  continuously  with  time  although 
they  can  sustain  a  direction  and  speed  for  a  period  of  days. 

Except  for  rare  events,  the  transient  wind  conditions  computed 
with  the  model  cannot  usually  be  found  in  field  current  meter 
measurements . 

For  an  estuary  as  large  as  the  Potomac,  the  spatial  dis¬ 
tribution  of  wind  stress  may  also  be  important.  Winds  confined 
to  the  lower  30  nm  of  the  estuary  would  affect  circulation  up 
estuary  differently  than  if  along  the  whole  estuary.  For  the 
Potomac  there  is  some  evidence  that  there  is  little  difference 
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Figure  23.  Velocity  with  Wind  Minus  Velocity 
without  Wind  Averaged  for  Successive 
Tidal  Periods  Following  Initiation  of 
1.5  m/s  Downestuary  Wind 
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between  wind  conditions  between  Patuxent  Naval  Air  Station  and 
Quantico,  which  are  near  nm  19  and  60,  respectively  (Elliott,  1976). 
For  winds  down  the  estuary,  however,  a  three-hour  lag  was  found 
from  the  upper  to  the  lower  station,  indicating  that  surface 
stress  increased  progressively  with  time  and  distance  along  the 
estuary  rather  than  reaching  its  full  magnitude  at  once. 

Proper  representation  of  the  temporal  and  spatial  variation 
of  wind  is  a  problem  in  input  data  development  and  formulation  of 
surface  sheer  stress  rather  than  a  problem  in  the  overall  nu¬ 
merical  hydrodynamics.  The  latter,  however,  cannot  produce  com¬ 
plete  answers  until  the  wind  problem  is  properly  formulated. 

The  surface  shear  stress  relationships  as  presently  used  are 
directly  and  instantaneously  related  to  wind  speed  with  no  time 
lags  or  surface  wind  wave  build  up. 
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Table  3 

Longitudinal  Salinity  Profile  after  10-day  Model  Simulation  to  Stationary  State 


i 


Model  segments  are  delineated  in  Figure 
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